The influence of non-driving related tasks on take-over quality in highly automated driving - a multilevel meta-analysis
Environment Preparation and Data Import
Import of relevant Libraries.
library(readr)
library(tidyverse)
library(readxl)
library(meta)
library(metafor)
library(esc)
library(dplyr)
library(meta)
library(metafor)
library(metaviz)
library(onehot)
library(mltools)
library(caret)
library(plyr)
library(ggplot2)
library(devtools)
library(dmetar)
library(car)
library(rmarkdown)
library(writexl)
library(rmdformats)
library(extrafont)
library(forestplot)
library(rmeta)
library(jtools)
library(PerformanceAnalytics)
library(psych)
library(GGally)
font_install("fontcm")
if (!require("devtools")) {
install.packages("devtools")
}
devtools::install_github("MathiasHarrer/dmetar")
remotes::install_github("crsh/papaja")
options(scipen=999)Generate functions to be used throughout the analysis.
# Function to automatically calculate effect sizes
effectsizecalculator <- function (i) {
esc_mean_sd(grp1m = data$minttc[effectsizes$id_ctrl[i] == data$id],
grp1sd = data$minttc_sd[effectsizes$id_ctrl[i] == data$id],
grp1n = data$sample_size[effectsizes$id_ctrl[i] == data$id],
grp2m = data$minttc[effectsizes$id_exp[i] == data$id],
grp2sd = data$minttc_sd[effectsizes$id_exp[i] == data$id],
grp2n = data$sample_size[effectsizes$id_exp[i] == data$id],
es.type = 'g',
study = studies$authors_meta[effectsizes$study_id[i] == studies$id])
}
# Function to locate authors to studies
authors <- function (a) {
studies$authors[effectsizes$study_id[a] == studies$id]
}
# Function to calculate variance distribution
mlm.variance.distribution = function(x){
m = x
# Check class
if (!(class(m)[1] %in% c("rma.mv", "rma"))){
stop("x must be of class 'rma.mv'.")
}
# Check for three level model
if (m$sigma2s != 2){
stop("The model you provided does not seem to be a three-level model. This function can only be used for three-level models.")
}
# Get variance diagonal and calculate total variance
n = m$k.eff
vector.inv.var = 1/(diag(m$V))
sum.inv.var = sum(vector.inv.var)
sum.sq.inv.var = (sum.inv.var)^2
vector.inv.var.sq = 1/(diag(m$V)^2)
sum.inv.var.sq = sum(vector.inv.var.sq)
num = (n-1)*sum.inv.var
den = sum.sq.inv.var - sum.inv.var.sq
est.samp.var = num/den
# Calculate variance proportions
level1=((est.samp.var)/(m$sigma2[1]+m$sigma2[2]+est.samp.var)*100)
level2=((m$sigma2[2])/(m$sigma2[1]+m$sigma2[2]+est.samp.var)*100)
level3=((m$sigma2[1])/(m$sigma2[1]+m$sigma2[2]+est.samp.var)*100)
# Prepare df for return
Level=c("Level 1", "Level 2", "Level 3")
Variance=c(level1, level2, level3)
df.res=data.frame(Variance)
colnames(df.res) = c("% of total variance")
rownames(df.res) = Level
I2 = c("---", round(Variance[2:3], 2))
df.res = as.data.frame(cbind(df.res, I2))
totalI2 = Variance[2] + Variance[3]
# Generate plot
df1 = data.frame("Level" = c("Sampling Error", "Total Heterogeneity"),
"Variance" = c(df.res[1,1], df.res[2,1]+df.res[3,1]),
"Type" = rep(1,2))
df2 = data.frame("Level" = rownames(df.res),
"Variance" = df.res[,1],
"Type" = rep(2,3))
df = as.data.frame(rbind(df1, df2))
g = ggplot(df, aes(fill=Level, y=Variance, x=as.factor(Type))) +
coord_cartesian(ylim = c(0,1), clip = "off") +
geom_bar(stat="identity", position="fill", width = 1, color="black") +
scale_y_continuous(labels = scales::percent)+
theme(axis.title.x=element_blank(),
axis.text.y = element_text(color="black"),
axis.line.y = element_blank(),
axis.title.y=element_blank(),
axis.line.x = element_blank(),
axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.y = element_line(lineend = "round"),
legend.position = "none",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
legend.background = element_rect(linetype="solid",
colour ="black"),
legend.title = element_blank(),
legend.key.size = unit(0.75,"cm"),
axis.ticks.length=unit(.25, "cm"),
plot.margin = unit(c(1,3,1,1), "lines")) +
scale_fill_manual(values = c("darkseagreen3", "deepskyblue3", "darkseagreen2",
"deepskyblue1", "deepskyblue2")) +
# Add Annotation
# Total Variance
annotate("text", x = 1.5, y = 1.05,
label = paste("Total Variance:",
round(m$sigma2[1]+m$sigma2[2]+est.samp.var, 3))) +
# Sampling Error
annotate("text", x = 1, y = (df[1,2]/2+df[2,2])/100,
label = paste("Sampling Error Variance: \n", round(est.samp.var, 3)), size = 3) +
# Total I2
annotate("text", x = 1, y = ((df[2,2])/100)/2-0.02,
label = bquote("Total"~italic(I)^2*":"~.(round(df[2,2],2))*"%"), size = 3) +
annotate("text", x = 1, y = ((df[2,2])/100)/2+0.05,
label = paste("Variance not attributable \n to sampling error: \n", round(m$sigma2[1]+m$sigma2[2],3)), size = 3) +
# Level 1
annotate("text", x = 2, y = (df[1,2]/2+df[2,2])/100, label = paste("Level 1: \n",
round(df$Variance[3],2), "%", sep=""), size = 3) +
# Level 2
annotate("text", x = 2, y = (df[5,2]+(df[4,2]/2))/100,
label = bquote(italic(I)[Level2]^2*":"~.(round(df[4,2],2))*"%"), size = 3) +
# Level 3
annotate("text", x = 2, y = (df[5,2]/2)/100,
label = bquote(italic(I)[Level3]^2*":"~.(round(df[5,2],2))*"%"), size = 3)
print(df.res)
cat("Total I2: ", round(totalI2, 2), "% \n", sep="")
suppressWarnings(print(g))
invisible(df.res)
}
# Automatically spot outliers
spot.outliers.random<-function(data){
data<-data
Author<-data$studlab
lowerci<-data$lower
upperci<-data$upper
m.outliers<-data.frame(Author,lowerci,upperci)
te.lower<-data$lower.random
te.upper<-data$upper.random
dplyr::filter(m.outliers,upperci < te.lower)
dplyr::filter(m.outliers,lowerci > te.upper)
}
# Plot hat values
hat.plot <-function(fit) {
p <- length(coefficients(fit))
n <- length(fitted(fit))
plot(hatvalues(fit), main="Index Plot of Hat Values")
abline(h=c(2,3)*p/n, col="red", lty=2)
identify(1:n, hatvalues(fit), names(hatvalues(fit)))
}Import files.
studies <- read_excel('~/Documents/1TU Dresden/4. Semester/Master-Thesis/Automatisiertes Fahren - Metaanalyse/00 Arbeitsmaterialien/Daten_final/studies.xlsx')
data <- read_excel('~/Documents/1TU Dresden/4. Semester/Master-Thesis/Automatisiertes Fahren - Metaanalyse/00 Arbeitsmaterialien/Daten_final/dataset.xlsx')
data_num <- read_excel('~/Documents/1TU Dresden/4. Semester/Master-Thesis/Automatisiertes Fahren - Metaanalyse/00 Arbeitsmaterialien/Daten_final/dataset_num.xlsx')
effectsizes <- read_excel('~/Documents/1TU Dresden/4. Semester/Master-Thesis/Automatisiertes Fahren - Metaanalyse/00 Arbeitsmaterialien/Daten_final/effsize_calc.xlsx')
head(studies)## # A tibble: 6 x 7
## id authors authors_meta year title type Source
## <dbl> <chr> <chr> <dbl> <chr> <chr> <chr>
## 1 1 Feldhütter,… Feldhütter et… 2017 How the Durati… confe… Advances in Er…
## 2 2 Körber, Gol… Körber et al.… 2016 The influence … journ… Transportation…
## 3 3 Gold, Körbe… Gold et al., … 2016 Taking Over Co… journ… Human Factors:…
## 4 4 Wandtner, S… Wandtner et a… 2018 Effects of Non… journ… Human Factors:…
## 5 6 Radlmayr, F… Radlmayr et a… 2019 The Influence … confe… Advances in In…
## 6 7 Wan & Wu Wan et al., 2… 2018 The Effects of… journ… International …## # A tibble: 6 x 27
## id studynr effectsize_id exp_cond exp_cond_descri… ndt_description
## <dbl> <dbl> <chr> <dbl> <chr> <chr>
## 1 1 1 1 1 SuRT visual, manual
## 2 2 1 1 2 No NDT -
## 3 3 2 2 1 Young, No NDT, … -
## 4 4 2 3 2 Young, No NDT, … -
## 5 5 2 4 3 Young, No NDT, … -
## 6 6 2 2 4 Young, NDT, 0 v… acoustic, cogn…
## # … with 21 more variables: subject_group <dbl>, sample_size <dbl>,
## # minttc <dbl>, minttc_sd <dbl>, mtot <dbl>, mtot_sd <dbl>, age <dbl>,
## # age_sd <dbl>, sim <chr>, lad <dbl>, ndt_v <chr>, ndt_a <chr>, ndt_m <chr>,
## # ndt_c <chr>, hand <chr>, ndt_p <chr>, tor_p <chr>, dre <chr>, iru <chr>,
## # urg <chr>, tbtc <dbl>## # A tibble: 6 x 27
## id studynr effectsize_id exp_cond exp_cond_descri… ndt_description
## <dbl> <dbl> <chr> <dbl> <chr> <chr>
## 1 1 1 1 1 SuRT visual, manual
## 2 2 1 1 2 No NDT -
## 3 3 2 2 1 Young, No NDT, … -
## 4 4 2 3 2 Young, No NDT, … -
## 5 5 2 4 3 Young, No NDT, … -
## 6 6 2 2 4 Young, NDT, 0 v… acoustic, cogn…
## # … with 21 more variables: subject_group <dbl>, sample_size <dbl>,
## # minttc <dbl>, minttc_sd <dbl>, mtot <dbl>, mtot_sd <dbl>, age <dbl>,
## # age_sd <dbl>, sim <dbl>, lad <dbl>, ndt_v <dbl>, ndt_a <dbl>, ndt_m <dbl>,
## # ndt_c <dbl>, hand <dbl>, ndt_p <dbl>, tor_p <dbl>, dre <dbl>, iru <dbl>,
## # urg <dbl>, tbtc <dbl>## # A tibble: 6 x 10
## study_id effectsize_id id_exp id_ctrl effect_size standard_error variance
## <dbl> <dbl> <dbl> <dbl> <lgl> <lgl> <lgl>
## 1 1 1 1 2 NA NA NA
## 2 2 2 6 3 NA NA NA
## 3 2 3 7 4 NA NA NA
## 4 2 4 8 5 NA NA NA
## 5 2 5 12 9 NA NA NA
## 6 2 6 13 10 NA NA NA
## # … with 3 more variables: lower_CI <lgl>, upper_CI <lgl>, weight <lgl>Data Preprocessing
Change datatypes.
data$minttc <- sapply(data$minttc, as.numeric)
data$minttc_sd <- sapply(data$minttc_sd, as.numeric)
data$effectsize_id <- sapply(data$effectsize_id, as.numeric)Convert files to dataframe.
data <- as.data.frame(data)
studies <- as.data.frame(studies)
effectsizes <- as.data.frame(effectsizes)Convert moderators to factor.
data$sim <- factor(data$sim)
data$lad <-factor(data$lad)
data$ndt_v <- factor(data$ndt_v)
data$ndt_a <- factor(data$ndt_a)
data$ndt_m <- factor(data$ndt_m)
data$ndt_c <- factor(data$ndt_c)
data$hand <- factor(data$hand)
data$ndt_p <- factor(data$ndt_p)
data$tor_p <- factor(data$tor_p)
data$dre <- factor(data$dre)
data$iru <-factor(data$iru)
data$urg <- factor(data$urg)Effectsize Calculation
Calculate effect sizes for the studies.
Change effectsizes columns to numeric.
effectsizes$effect_size <- sapply(effectsizes$effect_size, as.numeric)
effectsizes$standard_error <- sapply(effectsizes$standard_error, as.numeric)
effectsizes$variance <- sapply(effectsizes$variance,as.numeric)
effectsizes$lower_CI <- sapply(effectsizes$lower_CI, as.numeric)
effectsizes$upper_CI <- sapply(effectsizes$upper_CI, as.numeric)
effectsizes$weight <- sapply(effectsizes$weight, as.numeric)Merge effectsizs with studies for details.
Merge effectsizes with data for moderators.
Remove unnecessary columns from dataframe.
drop <- c("effectsize_id.y","authors_meta","authors","id_exp","id_ctrl","info","es","studynr","exp_cond","k/e","subject_group","")
effectsizes = effectsizes[,!(names(effectsizes) %in% drop)]Rename column(s).
Show output of new effectsizes dataframe.
## 'data.frame': 35 obs. of 37 variables:
## $ study_id : num 1 2 2 2 2 2 2 3 3 3 ...
## $ effectsize_id : num 1 2 3 4 5 6 7 8 9 10 ...
## $ effect_size : num 0.76 0.269 0.571 0.516 1.127 ...
## $ standard_error : num 0.402 0.335 0.34 0.339 0.36 ...
## $ variance : num 0.161 0.112 0.116 0.115 0.13 ...
## $ lower_CI : num -0.027 -0.3875 -0.0959 -0.1486 0.4212 ...
## $ upper_CI : num 1.548 0.925 1.238 1.181 1.832 ...
## $ weight : num 6.2 8.92 8.63 8.7 7.72 ...
## $ totaln : num 27 36 36 36 36 36 36 71 71 71 ...
## $ measure : chr "g" "g" "g" "g" ...
## $ study : chr "Feldhütter et al., 2017" "Körber et al., 2016" "Körber et al., 2016" "Körber et al., 2016" ...
## $ year : num 2017 2016 2016 2016 2016 ...
## $ title : chr "How the Duration of Automated Driving Influences Take-Over Performance and Gaze Behavior" "The influence of age on the take-over of vehicle control in highly automated driving" "The influence of age on the take-over of vehicle control in highly automated driving" "The influence of age on the take-over of vehicle control in highly automated driving" ...
## $ type : chr "conference paper" "journal article" "journal article" "journal article" ...
## $ Source : chr "Advances in Ergonomic Design of Systems, Products and Processes" "Transportation Research Part F: Traffic Psychology and Behaviour" "Transportation Research Part F: Traffic Psychology and Behaviour" "Transportation Research Part F: Traffic Psychology and Behaviour" ...
## $ exp_cond_description: chr "SuRT" "Young, NDT, 0 vehicle/km" "Young, NDT, 10 vehicle/km" "Young, NDT, 20 vehicle/km" ...
## $ ndt_description : chr "visual, manual" "acoustic, cognitive" "acoustic, cognitive" "acoustic, cognitive" ...
## $ sample_size : num 12 18 18 18 18 18 18 35 35 35 ...
## $ minttc : num 2.31 2.46 1.56 1.3 2.61 ...
## $ minttc_sd : num 0.745 0.512 0.564 0.611 0.66 0.475 0.572 1.17 1.03 1.17 ...
## $ mtot : num 2.25 2.76 3.7 3.66 2.62 ...
## $ mtot_sd : num 0.855 0.88 0.97 1.24 1.29 1.41 1.1 1.09 1.22 1.14 ...
## $ age : num 24.2 23.3 23.3 23.3 66.7 ...
## $ age_sd : num 2.09 2.6 2.6 2.6 4.56 4.56 4.56 22.2 22.2 22.2 ...
## $ sim : Factor w/ 3 levels "high fidelity",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ lad : Factor w/ 2 levels "0","1": 2 2 2 2 2 2 2 2 2 2 ...
## $ ndt_v : Factor w/ 2 levels "no","yes": 2 1 1 1 1 1 1 1 1 1 ...
## $ ndt_a : Factor w/ 2 levels "no","yes": 1 2 2 2 2 2 2 2 2 2 ...
## $ ndt_m : Factor w/ 2 levels "no","yes": 2 1 1 1 1 1 1 1 1 1 ...
## $ ndt_c : Factor w/ 2 levels "no","yes": 1 2 2 2 2 2 2 2 2 2 ...
## $ hand : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
## $ ndt_p : Factor w/ 2 levels "no","yes": 2 2 2 2 2 2 2 2 2 2 ...
## $ tor_p : Factor w/ 1 level "yes": 1 1 1 1 1 1 1 1 1 1 ...
## $ dre : Factor w/ 2 levels "high","medium": 2 2 2 2 2 2 2 2 2 2 ...
## $ iru : Factor w/ 2 levels "no","yes": 1 1 2 2 1 2 2 1 2 2 ...
## $ urg : Factor w/ 3 levels "high","low","medium": 1 1 1 1 1 1 1 1 1 1 ...
## $ tbtc : num 6 7 7 7 7 7 7 7 7 7 ...print a copy of the effectsizes dataframe.
Descriptive Statistics
Descriptive statstics summary
## id studynr effectsize_id exp_cond
## Min. : 1.00 Min. : 1.000 Min. : 1.00 Min. : 1.000
## 1st Qu.: 14.75 1st Qu.: 2.750 1st Qu.: 7.00 1st Qu.: 2.000
## Median : 64.50 Median : 6.000 Median :45.00 Median : 4.000
## Mean : 52.29 Mean : 6.036 Mean :33.98 Mean : 4.234
## 3rd Qu.: 78.25 3rd Qu.: 8.000 3rd Qu.:58.00 3rd Qu.: 6.000
## Max. :110.00 Max. :12.000 Max. :74.00 Max. :12.000
## NA's :7 NA's :9
## exp_cond_description ndt_description subject_group sample_size
## Length:56 Length:56 Min. :1.000 Min. :12.00
## Class :character Class :character 1st Qu.:1.000 1st Qu.:18.00
## Mode :character Mode :character Median :1.000 Median :36.00
## Mean :1.464 Mean :32.54
## 3rd Qu.:2.000 3rd Qu.:39.00
## Max. :4.000 Max. :49.00
##
## minttc minttc_sd mtot mtot_sd
## Min. : 0.4375 Min. :0.0870 Min. :1.125 Min. :0.2400
## 1st Qu.: 1.5295 1st Qu.:0.5507 1st Qu.:1.833 1st Qu.:0.4250
## Median : 2.3315 Median :0.7345 Median :2.278 Median :0.8123
## Mean : 3.6273 Mean :1.3714 Mean :2.440 Mean :0.8025
## 3rd Qu.: 4.0499 3rd Qu.:1.1425 3rd Qu.:3.195 3rd Qu.:1.1200
## Max. :12.4960 Max. :7.7600 Max. :3.700 Max. :1.4400
## NA's :1 NA's :1
## age age_sd sim lad ndt_v ndt_a
## Min. :23.28 Min. : 2.090 high fidelity : 5 0:12 no :33 no :42
## 1st Qu.:25.73 1st Qu.: 4.560 low fidelity : 6 1:44 yes:23 yes:14
## Median :31.23 Median : 6.380 medium fidelity:45
## Mean :37.19 Mean : 7.527
## 3rd Qu.:45.00 3rd Qu.: 7.100
## Max. :71.18 Max. :22.200
## NA's :1 NA's :4
## ndt_m ndt_c hand ndt_p tor_p dre iru urg
## no :48 no :46 no :49 no :21 yes:56 high :13 no :22 high :42
## yes: 8 yes:10 yes: 7 yes:35 medium:43 yes:34 low : 8
## medium: 6
##
##
##
##
## tbtc
## Min. : 3.000
## 1st Qu.: 6.000
## Median : 7.000
## Mean : 8.357
## 3rd Qu.: 7.750
## Max. :20.000
## Histograms
Histogram of minTTC
ggplot(data=data, aes(minttc)) +
geom_histogram(breaks=seq(0, 15, by=0.4),
col="black",
fill="black",
alpha = .2,
binwidth = 5) +
labs(title="Histogram of minTTC") +
labs(x="minTTC (s)", y="frequency") +
theme_apa()
Histogram of minTTC standard deviation
ggplot(data=data, aes(minttc_sd)) +
geom_histogram(breaks=seq(0, 10, by=0.3),
col="black",
fill="black",
alpha = .2) +
labs(title="Histogram of minTTC standard deviance") +
labs(x="standard deviance of minTTC (s)", y="frequency") +
theme_apa()
Histogram of effect sizes
ggplot(data=effectsizes, aes(x=effect_size)) +
geom_histogram(breaks=seq(-4, 4, by=0.2),
col="black",
fill="black",
alpha = .2) +
labs(title="Histogram of calculated effect sizes") +
labs(x="effect size", y="frequency") +
theme_apa()
Histogram of mTOT
ggplot(data=data, aes(x=mtot)) +
geom_histogram(breaks=seq(0, 5, by=0.2),
col="black",
fill="black",
alpha = .2) +
labs(title="Histogram of meanTOT") +
labs(x="meanTOT (s)", y="frequency") +
theme_apa()## Warning: Removed 1 rows containing non-finite values (stat_bin).
Histogram of mTOT standard deviation
ggplot(data=data, aes(x=mtot_sd)) +
geom_histogram(breaks=seq(0, 5, by=0.2),
col="black",
fill="black",
alpha = .2) +
labs(title="Histogram of meanTOT standard deviance") +
labs(x="standard deviance of meanTOT (s)", y="frequency") +
theme_apa()## Warning: Removed 1 rows containing non-finite values (stat_bin).
Histogram of age
ggplot(data=data, aes(age)) +
geom_histogram(breaks=seq(0, 75, by=1),
col="black",
fill="black",
alpha = .2) +
geom_density(alpha=.2, fill="#FF6666") +
labs(title="Histogram of age") +
labs(x="age", y="frequency") +
theme_apa()
Correlation matrices
Pearson correlation
# prepare dataset
vcols <- data_num %>% select(minttc, ndt_v, ndt_a, ndt_m, ndt_c, hand, ndt_p, dre, iru, urg)
var.cor <- cor(vcols, method = "pearson")
print(var.cor)## minttc ndt_v ndt_a ndt_m ndt_c hand
## minttc 1.00000000 0.07000871 -0.30620395 -0.21934296 -0.25659431 -0.23232672
## ndt_v 0.07000871 1.00000000 -0.23052136 0.38528033 -0.38924947 0.34298103
## ndt_a -0.30620395 -0.23052136 1.00000000 -0.23570226 0.80757285 -0.09352195
## ndt_m -0.21934296 0.38528033 -0.23570226 1.00000000 -0.19034675 0.61721340
## ndt_c -0.25659431 -0.38924947 0.80757285 -0.19034675 1.00000000 -0.17622684
## hand -0.23232672 0.34298103 -0.09352195 0.61721340 -0.17622684 1.00000000
## ndt_p -0.16550810 0.64666979 0.44721360 0.31622777 0.36115756 0.29277002
## dre 0.03814629 0.31472198 -0.21977383 -0.10360238 -0.25636488 -0.07993097
## iru -0.38025321 -0.07166747 0.12666010 0.01492704 0.08864896 -0.13819750
## urg -0.72899178 -0.19873812 0.25630730 0.08054961 0.25298389 -0.09321806
## ndt_p dre iru urg
## minttc -0.16550810 0.03814629 -0.38025321 -0.72899178
## ndt_v 0.64666979 0.31472198 -0.07166747 -0.19873812
## ndt_a 0.44721360 -0.21977383 0.12666010 0.25630730
## ndt_m 0.31622777 -0.10360238 0.01492704 0.08054961
## ndt_c 0.36115756 -0.25636488 0.08864896 0.25298389
## hand 0.29277002 -0.07993097 -0.13819750 -0.09321806
## ndt_p 1.00000000 0.07644455 0.05664412 0.03820804
## dre 0.07644455 1.00000000 0.44229225 0.29833845
## iru 0.05664412 0.44229225 1.00000000 0.67452786
## urg 0.03820804 0.29833845 0.67452786 1.00000000Spearman correlation
## minttc ndt_v ndt_a ndt_m ndt_c hand
## minttc 1.000000000 -0.07410923 -0.22835664 -0.19101569 -0.21491000 -0.30065935
## ndt_v -0.074109234 1.00000000 -0.23052136 0.38528033 -0.38924947 0.34298103
## ndt_a -0.228356640 -0.23052136 1.00000000 -0.23570226 0.80757285 -0.09352195
## ndt_m -0.191015687 0.38528033 -0.23570226 1.00000000 -0.19034675 0.61721340
## ndt_c -0.214910002 -0.38924947 0.80757285 -0.19034675 1.00000000 -0.17622684
## hand -0.300659355 0.34298103 -0.09352195 0.61721340 -0.17622684 1.00000000
## ndt_p -0.252172591 0.64666979 0.44721360 0.31622777 0.36115756 0.29277002
## dre 0.007850451 0.31472198 -0.21977383 -0.10360238 -0.25636488 -0.07993097
## iru -0.322361626 -0.07166747 0.12666010 0.01492704 0.08864896 -0.13819750
## urg -0.298769282 -0.24598667 0.24917096 0.03333333 0.26648545 -0.20720736
## ndt_p dre iru urg
## minttc -0.252172591 0.007850451 -0.32236163 -0.298769282
## ndt_v 0.646669791 0.314721976 -0.07166747 -0.245986671
## ndt_a 0.447213595 -0.219773831 0.12666010 0.249170961
## ndt_m 0.316227766 -0.103602377 0.01492704 0.033333333
## ndt_c 0.361157559 -0.256364882 0.08864896 0.266485446
## hand 0.292770022 -0.079930969 -0.13819750 -0.207207356
## ndt_p 1.000000000 0.076444546 0.05664412 0.003011693
## dre 0.076444546 1.000000000 0.44229225 0.314260545
## iru 0.056644118 0.442292251 1.00000000 0.710526906
## urg 0.003011693 0.314260545 0.71052691 1.000000000Scatterplots
Scatterplot of minTTC and minTTC standard deviation
ggplot(data=data, aes(minttc, minttc_sd)) +
geom_point(alpha = 0.5,) +
xlab("minTTC (s)") +
ylab("minTTC (s) standard deviance") +
theme_apa()
Scatterplot of and mTOT and mTOT standard deviation
ggplot(data=data, aes(mtot, mtot_sd)) +
geom_point(alpha = 0.5,) +
xlab("mTOT (s)") +
ylab("mTOT (s) standard deviance") +
theme_apa()## Warning: Removed 1 rows containing missing values (geom_point).
Scatterplot of meanTOT and minTTC
ggplot(data=data, aes(minttc, mtot)) +
geom_point(alpha = 0.5,) +
xlab("minTTC (s)") +
ylab("meanTOT (s)") +
theme_apa()## Warning: Removed 1 rows containing missing values (geom_point).
Scatterplot of meanTOT and minTTC with a fitted least squares line
## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 1 rows containing non-finite values (stat_smooth).## Warning: Removed 1 rows containing missing values (geom_point).
Scatterplot of TBTB and minTTC
Scatterplot of TBTC and minTTC with a fitted least squares line
Scatterplot Matrix
Scattterplot Matrix for continuous variables
mcols <- data %>% select(minttc, minttc_sd, mtot, mtot_sd, tbtc)
#pairs(~minttc+minttc_sd+mtot+mtot_sd+tbtc,data=data, pch=1)
ggpairs(mcols, axisLabels="show",
upper = list(continuous = wrap("cor", method="spearman", size = 3)),
lower = list(continuous = "points"))+
theme_apa()
Group comparisons
Group Comparisons for the dependent variable minTTC
# Group comparison minTTC - presence of NDT
count_ndt_p <- data %>% count("ndt_p")
minttc_compar_ndt_p <- data %>% group_by(ndt_p) %>% summarise_at(vars(minttc, minttc_sd), list(mean))
freq_ndt_p <- merge(minttc_compar_ndt_p,count_ndt_p)
print(freq_ndt_p)## ndt_p minttc minttc_sd freq
## 1 no 4.307712 1.715878 21
## 2 yes 3.219092 1.164652 35# Group comparison minTTC - visual task
count_ndt_v <- data %>% count("ndt_v")
minttc_compar_ndt_v <- data %>% group_by(ndt_v) %>% summarise_at(vars(minttc, minttc_sd), list(mean))
freq_ndt_v <- merge(minttc_compar_ndt_v,count_ndt_v)
print(freq_ndt_v)## ndt_v minttc minttc_sd freq
## 1 no 3.441214 1.335237 33
## 2 yes 3.894353 1.423194 23# Group comparison minTTC - acoustic task
count_ndt_a <- data %>% count("ndt_a")
minttc_compar_ndt_a <- data %>% group_by(ndt_a) %>% summarise_at(vars(minttc, minttc_sd), list(mean))
freq_ndt_a <- merge(minttc_compar_ndt_a,count_ndt_a)
print(freq_ndt_a)## ndt_a minttc minttc_sd freq
## 1 no 4.190265 1.6022124 42
## 2 yes 1.938503 0.6788112 14# Group comparison minTTC - manual task
count_ndt_m <- data %>% count("ndt_m")
minttc_compar_ndt_m <- data %>% group_by(ndt_m) %>% summarise_at(vars(minttc, minttc_sd), list(mean))
freq_ndt_m <- merge(minttc_compar_ndt_m,count_ndt_m)
print(freq_ndt_m)## ndt_m minttc minttc_sd freq
## 1 no 3.912466 1.5079731 48
## 2 yes 1.916476 0.5516965 8# Group comparison minTTC - cognitive task
count_ndt_c <- data %>% count("ndt_c")
minttc_compar_ndt_c <- data %>% group_by(ndt_c) %>% summarise_at(vars(minttc, minttc_sd), list(mean))
freq_ndt_c <- merge(minttc_compar_ndt_c,count_ndt_c)
print(freq_ndt_c)## ndt_c minttc minttc_sd freq
## 1 no 4.008286 1.5084548 46
## 2 yes 1.874904 0.7407357 10# Group comparison minTTC - handheld task
count_hand <- data %>% count("hand")
minttc_compar_hand <- data %>% group_by(hand) %>% summarise_at(vars(minttc, minttc_sd), list(mean))
freq_hand <- merge(minttc_compar_hand,count_hand)
print(freq_hand)## hand minttc minttc_sd freq
## 1 no 3.906941 1.478230 49
## 2 yes 1.670010 0.623289 7# Group comparison minTTC - urgency(urg)
count_urg <- data %>% count("urg")
minttc_compar_urg <- data %>% group_by(urg) %>% summarise_at(vars(minttc, minttc_sd), list(mean))
freq_urg <- merge(minttc_compar_urg,count_urg)
print(freq_urg)## urg minttc minttc_sd freq
## 1 high 2.664980 0.699459 42
## 2 low 10.340125 5.342375 8
## 3 medium 1.413333 0.780000 6# Group comparison minTTC - driver response (dre)
count_dre <- data %>% count("dre")
minttc_compar_dre <- data %>% group_by(dre) %>% summarise_at(vars(minttc, minttc_sd), list(mean))
freq_dre <- merge(minttc_compar_dre,count_dre)
print(freq_dre)## dre minttc minttc_sd freq
## 1 high 3.848241 0.8312336 13
## 2 medium 3.560536 1.5346568 43# Group comparison minTTC - interaction with other road users (iru)
count_iru <- data %>% count("iru")
minttc_compar_iru <- data %>% group_by(iru) %>% summarise_at(vars(minttc, minttc_sd), list(mean))
freq_iru <- merge(minttc_compar_iru,count_iru)
print(freq_iru)## iru minttc minttc_sd freq
## 1 no 5.132591 2.4351364 22
## 2 yes 2.653329 0.6830376 34Group Comparisons for mTOT
count_ndt_p <- data %>% count("ndt_p")
mtot_compar_ndt_p <- data %>% group_by(ndt_p) %>% summarise_at(vars(mtot, mtot_sd), list(mean), na.rm=TRUE)
freq_ndt_p <- merge(mtot_compar_ndt_p,count_ndt_p)
print(freq_ndt_p)## ndt_p mtot mtot_sd freq
## 1 no 2.252643 0.8036511 21
## 2 yes 2.546419 0.8017740 35# Group comparison mTOT - visual task
count_ndt_v <- data %>% count("ndt_v")
mtot_compar_ndt_v <- data %>% group_by(ndt_v) %>% summarise_at(vars(mtot, mtot_sd), list(mean), na.rm=TRUE)
freq_ndt_v <- merge(mtot_compar_ndt_v,count_ndt_v)
print(freq_ndt_v)## ndt_v mtot mtot_sd freq
## 1 no 2.502813 0.8866890 33
## 2 yes 2.351632 0.6852637 23# Group comparison mTOT - acoustic task
count_ndt_a <- data %>% count("ndt_a")
mtot_compar_ndt_a <- data %>% group_by(ndt_a) %>% summarise_at(vars(mtot, mtot_sd), list(mean), na.rm=TRUE)
freq_ndt_a <- merge(mtot_compar_ndt_a,count_ndt_a)
print(freq_ndt_a)## ndt_a mtot mtot_sd freq
## 1 no 2.299061 0.7363836 42
## 2 yes 2.851147 0.9959558 14# Group comparison mTOT - manual task
count_ndt_m <- data %>% count("ndt_m")
mtot_compar_ndt_m <- data %>% group_by(ndt_m) %>% summarise_at(vars(mtot, mtot_sd), list(mean), na.rm=TRUE)
freq_ndt_m <- merge(mtot_compar_ndt_m,count_ndt_m)
print(freq_ndt_m)## ndt_m mtot mtot_sd freq
## 1 no 2.452222 0.7991150 48
## 2 yes 2.365390 0.8220884 8# Group comparison mTOT - cognitive task
count_ndt_c <- data %>% count("ndt_c")
mtot_compar_ndt_c <- data %>% group_by(ndt_c) %>% summarise_at(vars(mtot, mtot_sd), list(mean), na.rm=TRUE)
freq_ndt_c <- merge(mtot_compar_ndt_c,count_ndt_c)
print(freq_ndt_c)## ndt_c mtot mtot_sd freq
## 1 no 2.291669 0.7380329 46
## 2 yes 3.105242 1.0923630 10# Group comparison mTOT - handheld task
count_hand <- data %>% count("hand")
mtot_compar_hand <- data %>% group_by(hand) %>% summarise_at(vars(mtot, mtot_sd), list(mean), na.rm=TRUE)
freq_hand <- merge(mtot_compar_hand,count_hand)
print(freq_hand)## hand mtot mtot_sd freq
## 1 no 2.397528 0.7855185 49
## 2 yes 2.728025 0.9186033 7# Group comparison mTOT - urgency(urg)
count_urg <- data %>% count("urg")
mtot_compar_urg <- data %>% group_by(urg) %>% summarise_at(vars(mtot, mtot_sd), list(mean), na.rm=TRUE)
freq_urg <- merge(mtot_compar_urg,count_urg)
print(freq_urg)## urg mtot mtot_sd freq
## 1 high 2.442964 0.8294417 42
## 2 low 1.924500 0.4260000 8
## 3 medium 3.103333 1.1200000 6# Group comparison mTOT - driver response (dre)
count_dre <- data %>% count("dre")
mtot_compar_dre <- data %>% group_by(dre) %>% summarise_at(vars(mtot, mtot_sd), list(mean), na.rm=TRUE)
freq_dre <- merge(mtot_compar_dre,count_dre)
print(freq_dre)## dre mtot mtot_sd freq
## 1 high 1.742769 0.4050000 13
## 2 medium 2.655275 0.9254788 43# Group comparison mTOT - interaction with other road users (iru)
count_iru <- data %>% count("iru")
mtot_compar_iru <- data %>% group_by(iru) %>% summarise_at(vars(mtot, mtot_sd), list(mean), na.rm=TRUE)
freq_iru <- merge(mtot_compar_iru,count_iru)
print(freq_iru)## iru mtot mtot_sd freq
## 1 no 2.440727 0.8176364 22
## 2 yes 2.438834 0.7923367 34Mean difference in minTTC (s) between “NDT” and “no NDT” condition
means <- data %>% group_by(`ndt_p`) %>% summarise_at(vars(minttc, minttc_sd), list(mean))
means$upper <- means$minttc + (0.5 * means$minttc_sd)
means$lower <- means$minttc - (0.5 * means$minttc_sd)
compar_minttc <- data %>% group_by(`ndt_p`) %>% summarise_at(vars(minttc, minttc_sd), list(mean))
compar_minttc$upper <- compar_minttc$minttc + (0.5 * means$minttc_sd)
compar_minttc$lower <- compar_minttc$minttc - (0.5 * means$minttc_sd)
ggplot(compar_minttc, aes(x = `ndt_p`, y=`minttc`)) +
geom_col(width=0.6, fill=c("grey", "white"), col="black") +
geom_errorbar(data=compar_minttc, mapping=aes(x=`ndt_p`, ymin=upper, ymax=lower, width=0.1)) +
geom_text(aes(label = format(minttc, digits = 3), y = 0.7)) +
scale_x_discrete(labels=c("NDT", "No NDT")) +
xlab("Engaged vs not engaged in a NDT while driving highly automated") +
ylab("minTTC (s)") +
annotate("text", x = 1.5, y=5, label = "mean difference = 1.09") +
theme_apa()
Mean difference in mTOT (s) between “NDT” and “no NDT” condition
compar_mtot <- data %>% group_by(`ndt_p`) %>% summarise_at(vars(mtot, mtot_sd), list(mean), na.rm=TRUE)
compar_mtot$upper <- compar_mtot$mtot + (0.5 * compar_mtot$mtot_sd)
compar_mtot$lower <- compar_mtot$mtot - (0.5 * compar_mtot$mtot_sd)
ggplot(compar_mtot, aes(x = `ndt_p`, y=`mtot`)) +
geom_col(width=0.6, fill=c("grey", "white"), col="black") +
geom_errorbar(data=compar_mtot, mapping=aes(x=`ndt_p`, ymin=upper, ymax=lower, width=0.1)) +
geom_text(aes(label = format(mtot, digits = 3), y = 0.7)) +
scale_x_discrete(labels=c("NDT", "No NDT")) +
xlab("Engaged vs not engaged in a NDT while driving highly automated") +
ylab("mTOT (s)") +
annotate("text", x = 1.5, y=2.8, label = "mean difference = 0.3") +
theme_apa()
## id studynr effectsize_id exp_cond
## Min. : 1.00 Min. : 1.000 Min. : 1.00 Min. : 1.000
## 1st Qu.: 14.75 1st Qu.: 2.750 1st Qu.: 7.00 1st Qu.: 2.000
## Median : 64.50 Median : 6.000 Median :45.00 Median : 4.000
## Mean : 52.29 Mean : 6.036 Mean :33.98 Mean : 4.234
## 3rd Qu.: 78.25 3rd Qu.: 8.000 3rd Qu.:58.00 3rd Qu.: 6.000
## Max. :110.00 Max. :12.000 Max. :74.00 Max. :12.000
## NA's :7 NA's :9
## exp_cond_description ndt_description subject_group sample_size
## Length:56 Length:56 Min. :1.000 Min. :12.00
## Class :character Class :character 1st Qu.:1.000 1st Qu.:18.00
## Mode :character Mode :character Median :1.000 Median :36.00
## Mean :1.464 Mean :32.54
## 3rd Qu.:2.000 3rd Qu.:39.00
## Max. :4.000 Max. :49.00
##
## minttc minttc_sd mtot mtot_sd
## Min. : 0.4375 Min. :0.0870 Min. :1.125 Min. :0.2400
## 1st Qu.: 1.5295 1st Qu.:0.5507 1st Qu.:1.833 1st Qu.:0.4250
## Median : 2.3315 Median :0.7345 Median :2.278 Median :0.8123
## Mean : 3.6273 Mean :1.3714 Mean :2.440 Mean :0.8025
## 3rd Qu.: 4.0499 3rd Qu.:1.1425 3rd Qu.:3.195 3rd Qu.:1.1200
## Max. :12.4960 Max. :7.7600 Max. :3.700 Max. :1.4400
## NA's :1 NA's :1
## age age_sd sim lad ndt_v ndt_a
## Min. :23.28 Min. : 2.090 high fidelity : 5 0:12 no :33 no :42
## 1st Qu.:25.73 1st Qu.: 4.560 low fidelity : 6 1:44 yes:23 yes:14
## Median :31.23 Median : 6.380 medium fidelity:45
## Mean :37.19 Mean : 7.527
## 3rd Qu.:45.00 3rd Qu.: 7.100
## Max. :71.18 Max. :22.200
## NA's :1 NA's :4
## ndt_m ndt_c hand ndt_p tor_p dre iru urg
## no :48 no :46 no :49 no :21 yes:56 high :13 no :22 high :42
## yes: 8 yes:10 yes: 7 yes:35 medium:43 yes:34 low : 8
## medium: 6
##
##
##
##
## tbtc
## Min. : 3.000
## 1st Qu.: 6.000
## Median : 7.000
## Mean : 8.357
## 3rd Qu.: 7.750
## Max. :20.000
## Two-Level meta-analytic Model
Build a two-level model whose results then can be compared to the three-level model.
##
## Random-Effects Model (k = 35; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 0.689 (SE = 0.185)
## tau (square root of estimated tau^2 value): 0.830
## I^2 (total heterogeneity / total variability): 91.30%
## H^2 (total variability / sampling variability): 11.50
##
## Test for Heterogeneity:
## Q(df = 34) = 333.768, p-val < .001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.329 0.148 2.225 0.026 0.039 0.619 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## estimate ci.lb ci.ub
## tau^2 0.689 0.430 1.281
## tau 0.830 0.656 1.132
## I^2(%) 91.302 86.768 95.125
## H^2 11.497 7.557 20.514Three-Level meta-analytic Model
Intercept-Only Model
Estimate the overall effect by fitting an intercept-only model.
overall <- rma.mv(effect_size, variance, random = list(~ 1 | effectsize_id, ~ 1 | study_id), tdist=
TRUE, data=effectsizes)
# Request a print of the results stored in the object # ‘‘overall’’ in three digits.
summary(overall, digits=3)##
## Multivariate Meta-Analysis Model (k = 35; method: REML)
##
## logLik Deviance AIC BIC AICc
## -44.096 88.192 94.192 98.771 94.992
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.689 0.830 35 no effectsize_id
## sigma^2.2 0.000 0.000 9 no study_id
##
## Test for Heterogeneity:
## Q(df = 34) = 333.768, p-val < .001
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## 0.329 0.148 2.225 0.033 0.028 0.629 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The estimated effect is the same as for the two-level model, but the confidence intervals differ.
Variance Distribution
Determine how the total variance is distributed over the three levels of the meta-analytic model. Print the results in percentages on screen.
## [1] 8.697861## [1] 91.30214## [1] 0.0000001483077## % of total variance I2
## Level 1 8.6978606518637 ---
## Level 2 91.3021391998285 91.3
## Level 3 0.0000001483077 0## Total I2: 91.3%
Calculation of the Intraclass Correlation (ICC)
ICC describes how strongly units in the same group resemble each other.
## [1] 0ICC = 0 indicates, that the values within clusters (=studies) are not similar. Because of the underlying research question, this was to be expected.
Profile Likelihood Plot
σ_1^2 (variance at the effectsize level) and σ_2^2 (variance at the study level) were fixed at different values. For each value of σ_1^2 and σ_2^2, the likelihood over the remaining model parameters, such as the fixed effects, was estimated. This means it was estimated how likely the values of these parameters are given the observed data. Less negative values of the logarithm indicate a higher likelihood than more negative values.
##
|
| | 0%
|
|==== | 5%
|
|======= | 10%
|
|========== | 15%
|
|============== | 20%
|
|================== | 25%
|
|===================== | 30%
|
|======================== | 35%
|
|============================ | 40%
|
|================================ | 45%
|
|=================================== | 50%
|
|====================================== | 55%
|
|========================================== | 60%
|
|============================================== | 65%
|
|================================================= | 70%
|
|==================================================== | 75%
|
|======================================================== | 80%
|
|============================================================ | 85%
|
|=============================================================== | 90%
|
|================================================================== | 95%
|
|======================================================================| 100%##
|
| | 0%
|
|==== | 5%
|
|======= | 10%
|
|========== | 15%
|
|============== | 20%
|
|================== | 25%
|
|===================== | 30%
|
|======================== | 35%
|
|============================ | 40%
|
|================================ | 45%
|
|=================================== | 50%
|
|====================================== | 55%
|
|========================================== | 60%
|
|============================================== | 65%
|
|================================================= | 70%
|
|==================================================== | 75%
|
|======================================================== | 80%
|
|============================================================ | 85%
|
|=============================================================== | 90%
|
|================================================================== | 95%
|
|======================================================================| 100%
The likelihood are highest for the values of σ_1^2 and σ_2^2 that had been estimated in the original model. Since profiling for variance components is done for non-negative value of the variance, the line is highest at this point in the second plot. We can be “fairly confident” that our meta-analytic models could identify the variance components (Viechtbauer, 2017).
Confidence Intervals for Variance
##
## estimate ci.lb ci.ub
## sigma^2.1 0.6890 0.4095 1.2176
## sigma.1 0.8301 0.6399 1.1034
##
## estimate ci.lb ci.ub
## sigma^2.2 0.0000 0.0000 0.2900
## sigma.2 0.0000 0.0000 0.5385The confindence intervals are very wide. For 𝜎22 hence, while the most likely value is 0, other values above 0 cannot be rejected either.
Two-Level Model without within-study variance
Build a two-level model without within-study variance.
modelnovar2 <- rma.mv(effect_size, variance, random = list(~ 1 | effectsize_id, ~ 1 | study_id),
sigma2=c(0,NA), tdist=TRUE, data=effectsizes)
# Request a print of the results stored in the object # ‘‘movelnovar2’’ in three digits.
summary(modelnovar2, digits=3)##
## Multivariate Meta-Analysis Model (k = 35; method: REML)
##
## logLik Deviance AIC BIC AICc
## -144.225 288.450 292.450 295.503 292.837
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.000 0.000 35 yes effectsize_id
## sigma^2.2 0.059 0.244 9 no study_id
##
## Test for Heterogeneity:
## Q(df = 34) = 333.768, p-val < .001
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## 0.315 0.097 3.264 0.003 0.119 0.511 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Perform a likelihood-ratio-test to determine the significance of the within-study (level 2) variance.
##
## df AIC BIC AICc logLik LRT pval QE
## Full 3 94.1919 98.7710 94.9919 -44.0960 333.7676
## Reduced 2 292.4499 295.5026 292.8370 -144.2249 200.2579 <.0001 333.7676The model does not perform better than the three-level model.
Two-level model without between-study variance
Build a two-level model without between-study variance.
modelnovar3 <- rma.mv(effect_size, variance, random = list(~ 1 | effectsize_id, ~ 1 | study_id),
sigma2=c(NA,0), tdist=TRUE, data=effectsizes)
# Request a print of the results stored in the object # ‘‘modelnovar3’’ in three digits.
summary(modelnovar3, digits=3)##
## Multivariate Meta-Analysis Model (k = 35; method: REML)
##
## logLik Deviance AIC BIC AICc
## -44.096 88.192 92.192 95.245 92.579
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.689 0.830 35 no effectsize_id
## sigma^2.2 0.000 0.000 9 yes study_id
##
## Test for Heterogeneity:
## Q(df = 34) = 333.768, p-val < .001
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## 0.329 0.148 2.225 0.033 0.028 0.629 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Perform a likelihood-ratio-test to determine the significance of the between-study (level 3) variance.
##
## df AIC BIC AICc logLik LRT pval QE
## Full 3 94.1919 98.7710 94.9919 -44.0960 333.7676
## Reduced 2 92.1919 95.2447 92.5790 -44.0960 0.0000 1.0000 333.7676This model does perform very similar to the three-level model. A third model for modeling differences between studies is not needed.
Moderator Analysis
NDT Modality
Test of different modalities independently
— acoustic NDT
ndt_a <- rma.mv(effect_size, variance, mods = ~ ndt_a, random = list(~1 | effectsize_id, ~ 1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes)
summary(ndt_a, digits=3)##
## Multivariate Meta-Analysis Model (k = 35; method: REML)
##
## logLik Deviance AIC BIC AICc
## -43.300 86.599 92.599 97.089 93.427
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.713 0.844 35 no effectsize_id
## sigma^2.2 0.000 0.000 9 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 33) = 333.302, p-val < .001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 33) = 0.002, p-val = 0.969
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.334 0.193 1.731 0.093 -0.058 0.726 .
## ndt_ayes -0.012 0.307 -0.039 0.969 -0.637 0.613
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1— manual NDT
ndt_m <- rma.mv(effect_size, variance, mods = ~ ndt_m, random = list(~1 | effectsize_id, ~ 1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes)
summary(ndt_m, digits=3)##
## Multivariate Meta-Analysis Model (k = 35; method: REML)
##
## logLik Deviance AIC BIC AICc
## -42.849 85.697 91.697 96.187 92.525
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.694 0.833 35 no effectsize_id
## sigma^2.2 0.000 0.000 9 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 33) = 331.804, p-val < .001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 33) = 0.901, p-val = 0.350
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.254 0.168 1.507 0.141 -0.089 0.596
## ndt_myes 0.338 0.357 0.949 0.350 -0.387 1.064
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1— cognitive NDT
ndt_c <- rma.mv(effect_size, variance, mods = ~ ndt_c, random = list(~1 | effectsize_id, ~ 1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes)
summary(ndt_c, digits=3)##
## Multivariate Meta-Analysis Model (k = 35; method: REML)
##
## logLik Deviance AIC BIC AICc
## -43.137 86.274 92.274 96.763 93.101
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.705 0.840 35 no effectsize_id
## sigma^2.2 0.000 0.000 9 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 33) = 330.616, p-val < .001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 33) = 0.317, p-val = 0.577
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.277 0.176 1.576 0.125 -0.081 0.635
## ndt_cyes 0.188 0.334 0.563 0.577 -0.491 0.866
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1— visual NDT
ndt_c <- rma.mv(effect_size, variance, mods = ~ ndt_v, random = list(~1 | effectsize_id, ~ 1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes)
summary(ndt_c, digits=3)##
## Multivariate Meta-Analysis Model (k = 35; method: REML)
##
## logLik Deviance AIC BIC AICc
## -43.151 86.303 92.303 96.792 93.130
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.707 0.841 35 no effectsize_id
## sigma^2.2 0.000 0.000 9 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 33) = 332.635, p-val < .001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 33) = 0.290, p-val = 0.594
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.215 0.259 0.830 0.413 -0.312 0.742
## ndt_vyes 0.171 0.317 0.539 0.594 -0.475 0.817
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1—handheld NDT
ndt_hand <- rma.mv(effect_size, variance, mods = ~ hand, random = list(~1 | effectsize_id, ~ 1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes)
summary(ndt_hand, digits=3)##
## Multivariate Meta-Analysis Model (k = 35; method: REML)
##
## logLik Deviance AIC BIC AICc
## -42.923 85.846 91.846 96.336 92.674
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.698 0.835 35 no effectsize_id
## sigma^2.2 0.000 0.000 9 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 33) = 331.826, p-val < .001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 33) = 0.767, p-val = 0.388
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.264 0.166 1.588 0.122 -0.074 0.602
## handyes 0.325 0.372 0.876 0.388 -0.431 1.082
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1All task modalities in one model.
ndt_modality <- rma.mv(effect_size, variance, mods = ~ ndt_a + ndt_m + ndt_c + ndt_v + hand, random = list(~1 | effectsize_id, ~ 1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes)
summary(ndt_modality, digits=3)##
## Multivariate Meta-Analysis Model (k = 35; method: REML)
##
## logLik Deviance AIC BIC AICc
## -36.163 72.326 86.326 95.897 91.659
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.626 0.791 35 no effectsize_id
## sigma^2.2 0.000 0.000 9 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 29) = 283.692, p-val < .001
##
## Test of Moderators (coefficients 2:6):
## F(df1 = 5, df2 = 29) = 1.705, p-val = 0.165
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt -1.479 0.709 -2.087 0.046 -2.929 -0.030 *
## ndt_ayes 0.135 0.506 0.266 0.792 -0.900 1.169
## ndt_myes 0.417 0.458 0.910 0.370 -0.520 1.354
## ndt_cyes 1.809 0.721 2.509 0.018 0.335 3.283 *
## ndt_vyes 1.638 0.669 2.449 0.021 0.270 3.006 *
## handyes 0.333 0.447 0.745 0.462 -0.581 1.248
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Urgency
One-hot encoding of urgency.
effectsizes <- effectsizes %>%
mutate(urg_low = ifelse(urg=='low', 1, 0),
urg_medium = ifelse(urg=='medium', 1, 0),
urg_high = ifelse(urg=='high', 1, 0))Determine the potential moderating effect of urgency as categorical moderator. Low urgency is chosen as the reference category.
urgency <- rma.mv(effect_size, variance, mods = ~ urg_medium + urg_high, random = list(~1 | effectsize_id, ~ 1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes)
summary(urgency, digits=3)##
## Multivariate Meta-Analysis Model (k = 35; method: REML)
##
## logLik Deviance AIC BIC AICc
## -42.364 84.727 92.727 98.590 94.209
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.731 0.855 35 no effectsize_id
## sigma^2.2 0.000 0.000 9 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 32) = 326.517, p-val < .001
##
## Test of Moderators (coefficients 2:3):
## F(df1 = 2, df2 = 32) = 0.152, p-val = 0.859
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.351 0.443 0.792 0.434 -0.552 1.253
## urg_medium 0.175 0.595 0.295 0.770 -1.037 1.387
## urg_high -0.064 0.477 -0.134 0.894 -1.036 0.908
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Determine the potential moderating effect of urgency as continuous moderator. High urgency is chosen as intercept category.
urgency_cont <- rma.mv(effect_size, variance, mods = ~ urg, random = list(~1 | effectsize_id, ~ 1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes)
summary(urgency_cont, digits=3)##
## Multivariate Meta-Analysis Model (k = 35; method: REML)
##
## logLik Deviance AIC BIC AICc
## -42.364 84.727 92.727 98.590 94.209
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.731 0.855 35 no effectsize_id
## sigma^2.2 0.000 0.000 9 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 32) = 326.517, p-val < .001
##
## Test of Moderators (coefficients 2:3):
## F(df1 = 2, df2 = 32) = 0.152, p-val = 0.859
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.287 0.177 1.621 0.115 -0.074 0.647
## urglow 0.064 0.477 0.134 0.894 -0.908 1.036
## urgmedium 0.239 0.435 0.550 0.586 -0.647 1.125
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Determine the potential moderating effect of time budget to collision.
tbtc <- rma.mv(effect_size, variance, mods = ~ tbtc, random = list(~ 1 | effectsize_id, ~1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes)
summary(tbtc, digits=3)##
## Multivariate Meta-Analysis Model (k = 35; method: REML)
##
## logLik Deviance AIC BIC AICc
## -43.264 86.528 92.528 97.017 93.355
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.711 0.843 35 no effectsize_id
## sigma^2.2 0.000 0.000 9 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 33) = 331.490, p-val < .001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 33) = 0.106, p-val = 0.747
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.247 0.293 0.842 0.406 -0.350 0.844
## tbtc 0.010 0.031 0.325 0.747 -0.054 0.074
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Exploratory Moderator Analyses
Simulator Fidelity
sim_fidelity <- rma.mv(effect_size, variance, mods = ~ sim, random = list(~1 | effectsize_id, ~ 1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes)
summary(sim_fidelity, digits=3)##
## Multivariate Meta-Analysis Model (k = 35; method: REML)
##
## logLik Deviance AIC BIC AICc
## -41.697 83.394 91.394 97.257 92.876
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.702 0.838 35 no effectsize_id
## sigma^2.2 0.000 0.000 9 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 32) = 322.433, p-val < .001
##
## Test of Moderators (coefficients 2:3):
## F(df1 = 2, df2 = 32) = 0.813, p-val = 0.453
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.775 0.447 1.735 0.092 -0.135 1.685 .
## simlow fidelity -0.249 0.593 -0.420 0.677 -1.457 0.959
## simmedium fidelity -0.552 0.479 -1.152 0.258 -1.528 0.424
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Driver response complexity
dre <- rma.mv(effect_size, variance, mods = ~ dre, random = list(~1 | effectsize_id, ~ 1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes)
summary(dre, digits=3)##
## Multivariate Meta-Analysis Model (k = 35; method: REML)
##
## logLik Deviance AIC BIC AICc
## -42.306 84.613 90.613 95.102 91.440
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.666 0.816 35 no effectsize_id
## sigma^2.2 0.000 0.000 9 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 33) = 315.691, p-val < .001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 33) = 2.079, p-val = 0.159
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt -0.019 0.282 -0.067 0.947 -0.592 0.554
## dremedium 0.474 0.329 1.442 0.159 -0.195 1.144
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Interaction with other road users
iru <- rma.mv(effect_size, variance, mods = ~ iru, random = list(~1 | effectsize_id, ~ 1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes)
summary(iru, digits=3)##
## Multivariate Meta-Analysis Model (k = 35; method: REML)
##
## logLik Deviance AIC BIC AICc
## -42.933 85.866 91.866 96.356 92.694
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.694 0.833 35 no effectsize_id
## sigma^2.2 0.000 0.000 9 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 33) = 323.312, p-val < .001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 33) = 0.751, p-val = 0.392
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.496 0.243 2.041 0.049 0.002 0.990 *
## iruyes -0.266 0.307 -0.867 0.392 -0.890 0.358
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Level of Automated Driving
lad <- rma.mv(effect_size, variance, mods = ~ lad, random = list(~1 | effectsize_id, ~ 1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes)
summary(lad, digits=3)##
## Multivariate Meta-Analysis Model (k = 35; method: REML)
##
## logLik Deviance AIC BIC AICc
## -43.316 86.632 92.632 97.121 93.459
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.714 0.845 35 no effectsize_id
## sigma^2.2 0.000 0.000 9 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 33) = 333.557, p-val < .001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 33) = 0.001, p-val = 0.980
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.322 0.330 0.975 0.337 -0.350 0.994
## lad1 0.009 0.371 0.025 0.980 -0.745 0.764
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Forest Plot
Forest Plot for individual effect sizes.
#sort effectsizes
effectsizes_sorted <- effectsizes[with(effectsizes, order(year, study, effect_size)),]
effectsize_num <- c(1:35)
# Forest Plot sorted
forest_sorted <- viz_forest(x = effectsizes_sorted[, c("effect_size", "standard_error")],
variant = "classic",
type = "standard",
study_table = effectsizes_sorted$study,
method = "REML",
annotate_CI = T,
table_headers = c("Study"),
xlab = "Hedges''g",
summary_table="Overall effect",
col = "darkblue",
summary_col = "firebrick",
table_layout = matrix(c(1, 2, 2, 3), nrow = 1))
print(forest_sorted)
Publication Bias
Visualization: Funnel Plot
Simple funnel plot.
Contour-enhanced funnel plot.
Contour-enhanced sunset funnel plot estimating the power of studies.
Trim-and-fill method
This is a non-parametric approach, where pseudo-studies are being imputed until the funnel plot symmetry is restored.
##
## Estimated number of missing studies on the left side: 10 (SE = 3.8535)
##
## Random-Effects Model (k = 45; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -64.7588 129.5176 133.5176 137.0860 133.8103
##
## tau^2 (estimated amount of total heterogeneity): 1.0060 (SE = 0.2313)
## tau (square root of estimated tau^2 value): 1.0030
## I^2 (total heterogeneity / total variability): 93.71%
## H^2 (total variability / sampling variability): 15.90
##
## Test for Heterogeneity:
## Q(df = 44) = 565.7152, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0060 0.1553 0.0384 0.9694 -0.2985 0.3104
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Testing the results of the the trim-and-fill method for the three-level model against the two-level model
viz_funnel(overall,
method="REML",
contours_col="Greys",
xlab="Hedges' g",
trim_and_fill = TRUE,
trim_and_fill_side = "left",
egger = FALSE)
No differences in the results can be detected.
Egger Regression Symmetry Test
Test for asymmetry of the funnel plot. It tests for the Y intercept = 0 from a linear regression of normalized effect estimate.
egger <- rma.mv(effect_size, variance, mods = standard_error, random = list(~ 1 | study_id, ~ 1 | effectsize_id), tdist= TRUE, data=effectsizes, method="REML")
summary(egger , digits=3)##
## Multivariate Meta-Analysis Model (k = 35; method: REML)
##
## logLik Deviance AIC BIC AICc
## -42.372 84.744 92.744 98.730 94.172
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.000 0.000 9 no study_id
## sigma^2.2 0.674 0.821 35 no effectsize_id
##
## Test for Residual Heterogeneity:
## QE(df = 33) = 325.381, p-val < .001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 33) = 1.872, p-val = 0.180
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt -0.597 0.692 -0.862 0.395 -2.005 0.811
## mods 3.433 2.509 1.368 0.180 -1.672 8.537
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The p-value is > 0.05, which indicates non-existence of publication bias.
Rank-Correlation Test
Correlates the standardized treatment effect with the variance of the treatment effect using Kendall’s tau as the measure of association. Since it is not available for three-level meta-analysis, it was calculated with the two-level model.
##
## Rank Correlation Test for Funnel Plot Asymmetry
##
## Kendall's tau = 0.1731, p = 0.1484The p-value is > 0.05, which indicates non-existence of publication bias.
Outlier & influence analysis
Leave-one-out analysis
Shows, how much Hedge’s g change after removing the studies individually.
leaveoneout <- viz_forest(x = effectsizes_sorted[, c("effect_size", "standard_error", "weight")],
variant = "classic",
type = "sensitivity",
study_table = effectsizes_sorted$study,
method = "REML",
annotate_CI = T,
table_headers = c("Study"),
xlab = "Hedges' g",
summary_table="Overall effect",
col = "darkblue",
summary_col = "firebrick",
table_layout = matrix(c(1, 2, 2, 3), nrow = 1))
print(leaveoneout)
Cook’s Distances
Result for each effect size or study shows how much the regression model would change when the value would be removed.
Cook’s Distances for each observed Outcome.
cook_ind <- cooks.distance(overall)
plot(cook_ind, type="o", pch=19, xlab="Effect sizes", ylab="Cook's Distance")
Cook’s Distances for each Cluster.
cook_cluster <- cooks.distance(overall, cluster=effectsizes$study_id)
plot(cook_cluster, type="o", pch=19, xlab="Studies", ylab="Cook's Distance", xaxt="n")
axis(side=1, at=seq_along(cook_cluster), labels=as.numeric(names(cook_cluster)))
Hat values
The hat values are the fitted values, the predictions made by the model for each observation.
## 1 2 3 4 5 6 7
## 0.02570633 0.02728570 0.02716063 0.02719011 0.02670611 0.02703940 0.02708251
## 8 9 10 11 12 13 14
## 0.02929797 0.02927713 0.02931434 0.02810648 0.02886847 0.02839469 0.02596490
## 15 16 17 18 19 20 21
## 0.02717960 0.02718419 0.02697977 0.02922146 0.02930504 0.02928871 0.02933955
## 22 23 24 25 26 27 28
## 0.02922106 0.02934942 0.02893951 0.02959745 0.02980240 0.02991593 0.02991766
## 29 30 31 32 33 34 35
## 0.02934960 0.02952701 0.02936674 0.02941902 0.02987211 0.02989404 0.02993496
Model without Outliers on effectsize level
Two-Level Model without Outliers on effectsize level
Removal of Outliers
Build a two-level model whose robustness then can be compared to the three-level model.
##
## Random-Effects Model (k = 30; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 0.083 (SE = 0.040)
## tau (square root of estimated tau^2 value): 0.288
## I^2 (total heterogeneity / total variability): 56.02%
## H^2 (total variability / sampling variability): 2.27
##
## Test for Heterogeneity:
## Q(df = 29) = 66.313, p-val < .001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.391 0.071 5.471 <.001 0.251 0.531 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## estimate ci.lb ci.ub
## tau^2 0.083 0.029 0.223
## tau 0.288 0.171 0.472
## I^2(%) 56.024 31.098 77.438
## H^2 2.274 1.451 4.432Three-Level Model without Outliers
Intercept-Only Model
Estimate the overall effect by fitting an intercept-only model.
Request a print of the results stored in the object ‘‘overall’’ in three digits.
##
## Multivariate Meta-Analysis Model (k = 30; method: REML)
##
## logLik Deviance AIC BIC AICc
## -14.708 29.416 35.416 39.518 36.376
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.083 0.288 30 no effectsize_id
## sigma^2.2 0.000 0.000 9 no study_id
##
## Test for Heterogeneity:
## Q(df = 29) = 66.313, p-val < .001
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## 0.391 0.071 5.471 <.001 0.245 0.537 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Variance Distribution
Determine how the total variance is distributed over the three levels of the meta-analytic model. Print the results in percentages on screen.
## [1] 43.97589## [1] 56.02411## [1] 0.0000002139203## % of total variance I2
## Level 1 43.9758852554402 ---
## Level 2 56.0241145306396 56.02
## Level 3 0.0000002139203 0## Total I2: 56.02%
Profile Likelihood Plots
##
|
| | 0%
|
|==== | 5%
|
|======= | 10%
|
|========== | 15%
|
|============== | 20%
|
|================== | 25%
|
|===================== | 30%
|
|======================== | 35%
|
|============================ | 40%
|
|================================ | 45%
|
|=================================== | 50%
|
|====================================== | 55%
|
|========================================== | 60%
|
|============================================== | 65%
|
|================================================= | 70%
|
|==================================================== | 75%
|
|======================================================== | 80%
|
|============================================================ | 85%
|
|=============================================================== | 90%
|
|================================================================== | 95%
|
|======================================================================| 100%##
|
| | 0%
|
|==== | 5%
|
|======= | 10%
|
|========== | 15%
|
|============== | 20%
|
|================== | 25%
|
|===================== | 30%
|
|======================== | 35%
|
|============================ | 40%
|
|================================ | 45%
|
|=================================== | 50%
|
|====================================== | 55%
|
|========================================== | 60%
|
|============================================== | 65%
|
|================================================= | 70%
|
|==================================================== | 75%
|
|======================================================== | 80%
|
|============================================================ | 85%
|
|=============================================================== | 90%
|
|================================================================== | 95%
|
|======================================================================| 100%
##
## estimate ci.lb ci.ub
## sigma^2.1 0.0827 0.0221 0.2015
## sigma.1 0.2876 0.1487 0.4489
##
## estimate ci.lb ci.ub
## sigma^2.2 0.0000 0.0000 0.1478
## sigma.2 0.0000 0.0000 0.3845Calculation of Intraclass Coefficient
## [1] 0Total sum of heterogeneity
## [1] 0.083Residuals
qqnorm(residuals(overall_wo,type="pearson"),main="QQ plot: residuals")
qqline(residuals(overall_wo,type="pearson"),col="red")
Two-Level Model without within-study variance
Build a two-level model without within-study variance.
modelnovar2_wo <- rma.mv(effect_size, variance, random = list(~ 1 | effectsize_id, ~ 1 | study_id),
sigma2=c(0,NA), tdist=TRUE, data=effectsizes_filtered)Request a print of the results stored in the object # ‘‘movelnovar2_wo’’ in three digits.
##
## Multivariate Meta-Analysis Model (k = 30; method: REML)
##
## logLik Deviance AIC BIC AICc
## -20.150 40.301 44.301 47.035 44.762
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.000 0.000 30 yes effectsize_id
## sigma^2.2 0.032 0.180 9 no study_id
##
## Test for Heterogeneity:
## Q(df = 29) = 66.313, p-val < .001
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## 0.382 0.080 4.791 <.001 0.219 0.545 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Perform a likelihood-ratio-test to determine the significance of the within-study (level2) variance.
##
## df AIC BIC AICc logLik LRT pval QE
## Full 3 35.4161 39.5180 36.3761 -14.7081 66.3129
## Reduced 2 44.3008 47.0353 44.7623 -20.1504 10.8846 0.0010 66.3129Two-Level Model without between-study variance
Build a two-level model without between-study variance
modelnovar3_wo <- rma.mv(effect_size, variance, random = list(~ 1 | effectsize_id, ~ 1 | study_id),
sigma2=c(NA,0), tdist=TRUE, data=effectsizes_filtered)
# Request a print of the results stored in the object # ‘‘modelnovar3’’ in three digits.
summary(modelnovar3_wo, digits=3)##
## Multivariate Meta-Analysis Model (k = 30; method: REML)
##
## logLik Deviance AIC BIC AICc
## -14.708 29.416 33.416 36.151 33.878
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.083 0.288 30 no effectsize_id
## sigma^2.2 0.000 0.000 9 yes study_id
##
## Test for Heterogeneity:
## Q(df = 29) = 66.313, p-val < .001
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## 0.391 0.071 5.471 <.001 0.245 0.537 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Perform a likelihood-ratio-test to determine the significance of the between-study (level3) variance.
##
## df AIC BIC AICc logLik LRT pval QE
## Full 3 35.4161 39.5180 36.3761 -14.7081 66.3129
## Reduced 2 33.4161 36.1507 33.8777 -14.7081 0.0000 1.0000 66.3129Moderator Analysis
NDT Modality
Determine the potential moderating effect of NDT modality.
Test of different modalities independently.
Acoustic NDT
ndt_a_wo <- rma.mv(effect_size, variance, mods = ~ ndt_a, random = list(~1 | effectsize_id, ~ 1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes_filtered)
summary(ndt_a_wo, digits=3)##
## Multivariate Meta-Analysis Model (k = 30; method: REML)
##
## logLik Deviance AIC BIC AICc
## -14.404 28.808 34.808 38.805 35.808
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.085 0.292 30 no effectsize_id
## sigma^2.2 0.000 0.000 9 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 28) = 64.873, p-val < .001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 28) = 0.610, p-val = 0.441
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.345 0.093 3.714 <.001 0.155 0.536 ***
## ndt_ayes 0.115 0.147 0.781 0.441 -0.186 0.416
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Manual NDT
ndt_m_wo <- rma.mv(effect_size, variance, mods = ~ ndt_m, random = list(~1 | effectsize_id, ~ 1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes_filtered)
summary(ndt_m_wo, digits=3)##
## Multivariate Meta-Analysis Model (k = 30; method: REML)
##
## logLik Deviance AIC BIC AICc
## -14.097 28.194 34.194 38.191 35.194
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.081 0.285 30 no effectsize_id
## sigma^2.2 0.000 0.000 9 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 28) = 63.156, p-val < .001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 28) = 1.221, p-val = 0.279
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.432 0.080 5.377 <.001 0.267 0.597 ***
## ndt_myes -0.191 0.172 -1.105 0.279 -0.544 0.163
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Cognitive NDT
ndt_c_wo <- rma.mv(effect_size, variance, mods = ~ ndt_c, random = list(~1 | effectsize_id, ~ 1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes_filtered)
summary(ndt_c_wo, digits=3)##
## Multivariate Meta-Analysis Model (k = 30; method: REML)
##
## logLik Deviance AIC BIC AICc
## -14.551 29.101 35.101 39.098 36.101
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.087 0.295 30 no effectsize_id
## sigma^2.2 0.000 0.000 9 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 28) = 65.911, p-val < .001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 28) = 0.257, p-val = 0.616
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.368 0.086 4.277 <.001 0.192 0.544 ***
## ndt_cyes 0.081 0.160 0.507 0.616 -0.247 0.409
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Handheld NDT
ndt_hand_wo <- rma.mv(effect_size, variance, mods = ~ hand, random = list(~1 | effectsize_id, ~ 1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes_filtered)
summary(ndt_hand_wo, digits=3)##
## Multivariate Meta-Analysis Model (k = 30; method: REML)
##
## logLik Deviance AIC BIC AICc
## -13.979 27.958 33.958 37.955 34.958
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.082 0.286 30 no effectsize_id
## sigma^2.2 0.000 0.000 9 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 28) = 63.148, p-val < .001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 28) = 1.548, p-val = 0.224
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.437 0.080 5.447 <.001 0.272 0.601 ***
## handyes -0.217 0.174 -1.244 0.224 -0.573 0.140
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1All modalities in one model.
ndt_modality_wo <- rma.mv(effect_size, variance, mods = ~ ndt_a + ndt_m + ndt_c + hand, random = list(~1 | effectsize_id, ~ 1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes_filtered)
summary(ndt_modality_wo, digits=3)##
## Multivariate Meta-Analysis Model (k = 30; method: REML)
##
## logLik Deviance AIC BIC AICc
## -13.715 27.430 39.430 46.743 44.097
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.097 0.311 30 no effectsize_id
## sigma^2.2 0.000 0.000 9 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 25) = 61.580, p-val < .001
##
## Test of Moderators (coefficients 2:5):
## F(df1 = 4, df2 = 25) = 0.474, p-val = 0.754
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.422 0.123 3.425 0.002 0.168 0.675 **
## ndt_ayes 0.128 0.261 0.492 0.627 -0.409 0.666
## ndt_myes -0.061 0.236 -0.260 0.797 -0.547 0.424
## ndt_cyes -0.100 0.275 -0.365 0.718 -0.666 0.465
## handyes -0.185 0.228 -0.809 0.426 -0.654 0.285
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Urgency
One-hot encoding of hand-coded urgency
effectsizes_filtered <- effectsizes_filtered %>%
mutate(urg_low = ifelse(urg=='low', 1, 0),
urg_medium = ifelse(urg=='medium', 1, 0),
urg_high = ifelse(urg=='high', 1, 0))Determine the potential moderating effect of hand-coded urgency as categorical moderator. No urgency is chosen as the reference category
urgency_wo <- rma.mv(effect_size, variance, mods = ~ urg_medium + urg_high, random = list(~1 | effectsize_id, ~ 1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes_filtered)
summary(urgency_wo, digits=3)##
## Multivariate Meta-Analysis Model (k = 30; method: REML)
##
## logLik Deviance AIC BIC AICc
## -14.416 28.832 36.832 42.015 38.650
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.090 0.300 30 no effectsize_id
## sigma^2.2 0.000 0.000 9 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 27) = 64.017, p-val < .001
##
## Test of Moderators (coefficients 2:3):
## F(df1 = 2, df2 = 27) = 0.373, p-val = 0.692
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.348 0.190 1.837 0.077 -0.041 0.738 .
## urg_medium 0.177 0.256 0.690 0.496 -0.349 0.702
## urg_high 0.016 0.210 0.077 0.939 -0.414 0.446
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Determine the potential moderating effect of urgency as continuous moderator. High urgency is chosen as inttercept category
urgency_cont_wo <- rma.mv(effect_size, variance, mods = ~ urg, random = list(~1 | effectsize_id, ~ 1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes_filtered)
summary(urgency_cont_wo, digits=3)##
## Multivariate Meta-Analysis Model (k = 30; method: REML)
##
## logLik Deviance AIC BIC AICc
## -14.416 28.832 36.832 42.015 38.650
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.090 0.300 30 no effectsize_id
## sigma^2.2 0.000 0.000 9 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 27) = 64.017, p-val < .001
##
## Test of Moderators (coefficients 2:3):
## F(df1 = 2, df2 = 27) = 0.373, p-val = 0.692
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.365 0.089 4.082 <.001 0.181 0.548 ***
## urglow -0.016 0.210 -0.077 0.939 -0.446 0.414
## urgmedium 0.160 0.194 0.828 0.415 -0.237 0.558
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Determine the potential moderating effect of time budget to collision.
tbtc_wo <- rma.mv(effect_size, variance, mods = ~ tbtc, random = list(~ 1 | effectsize_id, ~1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes_filtered)
summary(tbtc_wo, digits=3)##
## Multivariate Meta-Analysis Model (k = 30; method: REML)
##
## logLik Deviance AIC BIC AICc
## -14.770 29.539 35.539 39.536 36.539
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.089 0.298 30 no effectsize_id
## sigma^2.2 0.000 0.000 9 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 28) = 66.238, p-val < .001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 28) = 0.005, p-val = 0.944
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.383 0.143 2.675 0.012 0.090 0.676 *
## tbtc 0.001 0.014 0.070 0.944 -0.028 0.030
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Further Moderator Analyses
Simulator fidelity
sim_fidelity_wo <- rma.mv(effect_size, variance, mods = ~ sim, random = list(~1 | effectsize_id, ~ 1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes_filtered)
summary(sim_fidelity_wo, digits=3)##
## Multivariate Meta-Analysis Model (k = 30; method: REML)
##
## logLik Deviance AIC BIC AICc
## -12.745 25.489 33.489 38.672 35.307
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.071 0.267 30 no effectsize_id
## sigma^2.2 0.000 0.000 9 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 27) = 57.327, p-val < .001
##
## Test of Moderators (coefficients 2:3):
## F(df1 = 2, df2 = 27) = 2.148, p-val = 0.136
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.840 0.270 3.113 0.004 0.286 1.393 **
## simlow fidelity -0.315 0.314 -1.003 0.325 -0.959 0.329
## simmedium fidelity -0.521 0.281 -1.853 0.075 -1.098 0.056 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#der_wo <- rma.mv(effect_size, standard_error, mods = ~ der, random = list(~1 | effectsize_id, ~ 1 | study_id), tdist=TRUE, data=effectsizes_filtered)
#summary(der_wo, digits=3)Interaction with other road users
iru_wo <- rma.mv(effect_size, standard_error, mods = ~ iru, random = list(~1 | effectsize_id, ~ 1 | study_id), tdist=TRUE, data=effectsizes_filtered)
summary(iru_wo, digits=3)##
## Multivariate Meta-Analysis Model (k = 30; method: REML)
##
## logLik Deviance AIC BIC AICc
## -15.639 31.278 39.278 44.606 41.017
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.000 0.000 30 no effectsize_id
## sigma^2.2 0.000 0.000 9 no study_id
##
## Test for Residual Heterogeneity:
## QE(df = 28) = 16.990, p-val = 0.949
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 28) = 0.704, p-val = 0.408
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.478 0.141 3.393 0.002 0.189 0.766 **
## iruyes -0.157 0.187 -0.839 0.408 -0.541 0.227
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Level of automated driving
lad_wo <- rma.mv(effect_size, standard_error, mods = ~ lad, random = list(~1 | effectsize_id, ~ 1 | study_id), tdist=TRUE, data=effectsizes_filtered)
summary(lad_wo, digits=3)##
## Multivariate Meta-Analysis Model (k = 30; method: REML)
##
## logLik Deviance AIC BIC AICc
## -15.931 31.862 39.862 45.191 41.601
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.000 0.000 30 no effectsize_id
## sigma^2.2 0.000 0.000 9 no study_id
##
## Test for Residual Heterogeneity:
## QE(df = 28) = 17.467, p-val = 0.939
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 28) = 0.227, p-val = 0.638
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.317 0.177 1.789 0.084 -0.046 0.680 .
## lad1 0.099 0.208 0.476 0.638 -0.327 0.526
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Publication Bias
Visualization: Funnel Plot
Trim-and-fill method
##
## Estimated number of missing studies on the left side: 4 (SE = 3.6206)
##
## Random-Effects Model (k = 34; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -21.0215 42.0430 46.0430 49.0360 46.4430
##
## tau^2 (estimated amount of total heterogeneity): 0.1233 (SE = 0.0481)
## tau (square root of estimated tau^2 value): 0.3512
## I^2 (total heterogeneity / total variability): 64.74%
## H^2 (total variability / sampling variability): 2.84
##
## Test for Heterogeneity:
## Q(df = 33) = 91.1965, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.3120 0.0763 4.0904 <.0001 0.1625 0.4615 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Egger Regression Symmetry Test (eggers.test not working, bc ‘overall’ is no ‘meta’-type. Regtest is carried out with the two-level model)
##
## Regression Test for Funnel Plot Asymmetry
##
## model: mixed-effects meta-regression model
## predictor: standard error
##
## test for funnel plot asymmetry: z = 0.7768, p = 0.4373The p-value is > 0.05, publication bias can be rejected.
Another method for calculating Egger’s regression test - same result though.
egger <- rma.mv(effect_size, variance, mods = standard_error, random = list(~ 1 | study_id, ~ 1 | effectsize_id), tdist= TRUE, data=effectsizes_filtered, method="REML")
summary(egger , digits=3)##
## Multivariate Meta-Analysis Model (k = 30; method: REML)
##
## logLik Deviance AIC BIC AICc
## -14.349 28.698 36.698 42.027 38.437
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.000 0.000 9 no study_id
## sigma^2.2 0.083 0.289 30 no effectsize_id
##
## Test for Residual Heterogeneity:
## QE(df = 28) = 64.368, p-val < .001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 28) = 0.603, p-val = 0.444
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.127 0.347 0.366 0.717 -0.584 0.838
## mods 1.017 1.309 0.777 0.444 -1.664 3.697
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Rank-Correlation Test
##
## Rank Correlation Test for Funnel Plot Asymmetry
##
## Kendall's tau = 0.2046, p = 0.1170Model without Lin et al., 2020
Two-Level Model without Lin et al., 2020
Removal of Outliers
Build a two-level model whose robustness then can be compared to the three-level model.
##
## Random-Effects Model (k = 29; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 0.422 (SE = 0.134)
## tau (square root of estimated tau^2 value): 0.649
## I^2 (total heterogeneity / total variability): 85.62%
## H^2 (total variability / sampling variability): 6.96
##
## Test for Heterogeneity:
## Q(df = 28) = 150.829, p-val < .001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.433 0.132 3.289 0.001 0.175 0.691 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## estimate ci.lb ci.ub
## tau^2 0.422 0.246 0.936
## tau 0.649 0.496 0.967
## I^2(%) 85.623 77.679 92.966
## H^2 6.956 4.480 14.217Three-Level-Model without Lin et al., 2020
Intercept-Only Model
Estimate the overall effect by fitting an intercept-only model.
Request a print of the results stored in the object ‘‘overall’’ in three digits.
##
## Multivariate Meta-Analysis Model (k = 29; method: REML)
##
## logLik Deviance AIC BIC AICc
## -31.095 62.191 68.191 72.187 69.191
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.422 0.649 29 no effectsize_id
## sigma^2.2 0.000 0.000 8 no study_id
##
## Test for Heterogeneity:
## Q(df = 28) = 150.829, p-val < .001
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## 0.433 0.132 3.289 0.003 0.163 0.703 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Variance Distribution
Determine how the total variance is distributed over the three levels of the meta-analytic model. Print the results in percentages on screen.
## [1] 14.37685## [1] 85.62315## [1] 0.00000005084847## % of total variance I2
## Level 1 14.37684599661965 ---
## Level 2 85.62315395253188 85.62
## Level 3 0.00000005084847 0## Total I2: 85.62%
Profile Likelihood Plots
##
|
| | 0%
|
|==== | 5%
|
|======= | 10%
|
|========== | 15%
|
|============== | 20%
|
|================== | 25%
|
|===================== | 30%
|
|======================== | 35%
|
|============================ | 40%
|
|================================ | 45%
|
|=================================== | 50%
|
|====================================== | 55%
|
|========================================== | 60%
|
|============================================== | 65%
|
|================================================= | 70%
|
|==================================================== | 75%
|
|======================================================== | 80%
|
|============================================================ | 85%
|
|=============================================================== | 90%
|
|================================================================== | 95%
|
|======================================================================| 100%##
|
| | 0%
|
|==== | 5%
|
|======= | 10%
|
|========== | 15%
|
|============== | 20%
|
|================== | 25%
|
|===================== | 30%
|
|======================== | 35%
|
|============================ | 40%
|
|================================ | 45%
|
|=================================== | 50%
|
|====================================== | 55%
|
|========================================== | 60%
|
|============================================== | 65%
|
|================================================= | 70%
|
|==================================================== | 75%
|
|======================================================== | 80%
|
|============================================================ | 85%
|
|=============================================================== | 90%
|
|================================================================== | 95%
|
|======================================================================| 100%
##
## estimate ci.lb ci.ub
## sigma^2.1 0.4216 0.2168 0.8389
## sigma.1 0.6493 0.4656 0.9159
##
## estimate ci.lb ci.ub
## sigma^2.2 0.0000 0.0000 0.1884
## sigma.2 0.0000 0.0000 0.4341Calculation of Intraclass Coefficient
## [1] 0Total sum of heterogeneity
## [1] 0.422Residuals
qqnorm(residuals(overall_w8,type="pearson"),main="QQ plot: residuals")
qqline(residuals(overall_w8,type="pearson"),col="red")
Two-Level Model without within-study variance
Build a two-level model without within-study variance.
modelnovar2_w8 <- rma.mv(effect_size, variance, random = list(~ 1 | effectsize_id, ~ 1 | study_id),
sigma2=c(0,NA), tdist=TRUE, data=effectsizes_filtered2)Request a print of the results stored in the object # ‘‘movelnovar2’’ in three digits.
##
## Multivariate Meta-Analysis Model (k = 29; method: REML)
##
## logLik Deviance AIC BIC AICc
## -64.941 129.882 133.882 136.547 134.362
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.000 0.000 29 yes effectsize_id
## sigma^2.2 0.013 0.114 8 no study_id
##
## Test for Heterogeneity:
## Q(df = 28) = 150.829, p-val < .001
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## 0.390 0.066 5.868 <.001 0.254 0.526 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Perform a likelihood-ratio-test to determine the significance of the within-study (level2) variance.
##
## df AIC BIC AICc logLik LRT pval QE
## Full 3 68.1908 72.1874 69.1908 -31.0954 150.8287
## Reduced 2 133.8823 136.5467 134.3623 -64.9412 67.6915 <.0001 150.8287Two-Level Model without between-study variance
Build a two-level model without between-study variance
modelnovar3_w8 <- rma.mv(effect_size, variance, random = list(~ 1 | effectsize_id, ~ 1 | study_id),
sigma2=c(NA,0), tdist=TRUE, data=effectsizes_filtered2)
# Request a print of the results stored in the object # ‘‘modelnovar3’’ in three digits.
summary(modelnovar3_w8, digits=3)##
## Multivariate Meta-Analysis Model (k = 29; method: REML)
##
## logLik Deviance AIC BIC AICc
## -31.095 62.191 66.191 68.855 66.671
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.422 0.649 29 no effectsize_id
## sigma^2.2 0.000 0.000 8 yes study_id
##
## Test for Heterogeneity:
## Q(df = 28) = 150.829, p-val < .001
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## 0.433 0.132 3.289 0.003 0.163 0.703 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Perform a likelihood-ratio-test to determine the significance of the between-study (level3) variance.
##
## df AIC BIC AICc logLik LRT pval QE
## Full 3 68.1908 72.1874 69.1908 -31.0954 150.8287
## Reduced 2 66.1908 68.8552 66.6708 -31.0954 0.0000 1.0000 150.8287Moderator Analysis
NDT Modality
Determine the potential moderating effect of ndt modality.
Test of different modalities independently
ndt_a_w8 <- rma.mv(effect_size, variance, mods = ~ ndt_a, random = list(~1 | effectsize_id, ~ 1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes_filtered2)
summary(ndt_a_wo, digits=3)##
## Multivariate Meta-Analysis Model (k = 30; method: REML)
##
## logLik Deviance AIC BIC AICc
## -14.404 28.808 34.808 38.805 35.808
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.085 0.292 30 no effectsize_id
## sigma^2.2 0.000 0.000 9 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 28) = 64.873, p-val < .001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 28) = 0.610, p-val = 0.441
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.345 0.093 3.714 <.001 0.155 0.536 ***
## ndt_ayes 0.115 0.147 0.781 0.441 -0.186 0.416
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1— manual
ndt_m_w8 <- rma.mv(effect_size, variance, mods = ~ ndt_m, random = list(~1 | effectsize_id, ~ 1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes_filtered2)
summary(ndt_m_w8, digits=3)##
## Multivariate Meta-Analysis Model (k = 29; method: REML)
##
## logLik Deviance AIC BIC AICc
## -30.283 60.565 66.565 70.453 67.609
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.437 0.661 29 no effectsize_id
## sigma^2.2 0.000 0.000 8 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 27) = 150.814, p-val < .001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 27) = 0.440, p-val = 0.513
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.380 0.156 2.430 0.022 0.059 0.700 *
## ndt_myes 0.200 0.301 0.664 0.513 -0.418 0.818
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1— cognitive
ndt_c_w8 <- rma.mv(effect_size, variance, mods = ~ ndt_c, random = list(~1 | effectsize_id, ~ 1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes_filtered2)
summary(ndt_c_w8, digits=3)##
## Multivariate Meta-Analysis Model (k = 29; method: REML)
##
## logLik Deviance AIC BIC AICc
## -30.491 60.981 66.981 70.869 68.025
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.442 0.665 29 no effectsize_id
## sigma^2.2 0.000 0.000 8 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 27) = 150.647, p-val < .001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 27) = 0.023, p-val = 0.879
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.419 0.165 2.546 0.017 0.081 0.757 *
## ndt_cyes 0.044 0.285 0.153 0.879 -0.541 0.628
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1—handheld
ndt_hand_w8 <- rma.mv(effect_size, variance, mods = ~ hand, random = list(~1 | effectsize_id, ~ 1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes_filtered2)
summary(ndt_hand_w8, digits=3)##
## Multivariate Meta-Analysis Model (k = 29; method: REML)
##
## logLik Deviance AIC BIC AICc
## -30.335 60.669 66.669 70.557 67.713
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.439 0.662 29 no effectsize_id
## sigma^2.2 0.000 0.000 8 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 27) = 150.826, p-val < .001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 27) = 0.359, p-val = 0.554
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.388 0.154 2.518 0.018 0.072 0.704 *
## handyes 0.187 0.312 0.599 0.554 -0.453 0.826
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1All modalities in one model
ndt_modality_w8 <- rma.mv(effect_size, variance, mods = ~ ndt_a + ndt_m + ndt_c + hand, random = list(~1 | effectsize_id, ~ 1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes_filtered2)
summary(ndt_modality_w8, digits=3)##
## Multivariate Meta-Analysis Model (k = 29; method: REML)
##
## logLik Deviance AIC BIC AICc
## -27.622 55.244 67.244 74.312 72.185
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.473 0.688 29 no effectsize_id
## sigma^2.2 0.000 0.000 8 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 24) = 145.148, p-val < .001
##
## Test of Moderators (coefficients 2:5):
## F(df1 = 4, df2 = 24) = 0.524, p-val = 0.719
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.462 0.280 1.650 0.112 -0.116 1.040
## ndt_ayes -0.525 0.460 -1.141 0.265 -1.473 0.424
## ndt_myes -0.006 0.431 -0.014 0.989 -0.895 0.884
## ndt_cyes 0.525 0.449 1.171 0.253 -0.401 1.452
## handyes 0.197 0.399 0.493 0.627 -0.627 1.021
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Urgency
One-hot encoding of urgency
effectsizes_filtered2 <- effectsizes_filtered2 %>%
mutate(urg_low = ifelse(urg=='low', 1, 0),
urg_medium = ifelse(urg=='medium', 1, 0),
urg_high = ifelse(urg=='high', 1, 0))Determine the potential moderating effect of urgency as categorical moderator. No urgency is chosen as the reference category
urgency_w8 <- rma.mv(effect_size, variance, mods = ~ urg_medium + urg_high, random = list(~1 | effectsize_id, ~ 1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes_filtered2)
summary(urgency_w8, digits=3)##
## Multivariate Meta-Analysis Model (k = 29; method: REML)
##
## logLik Deviance AIC BIC AICc
## -29.872 59.743 67.743 72.776 69.648
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.464 0.681 29 no effectsize_id
## sigma^2.2 0.000 0.000 8 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 26) = 148.978, p-val < .001
##
## Test of Moderators (coefficients 2:3):
## F(df1 = 2, df2 = 26) = 0.068, p-val = 0.934
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.350 0.360 0.974 0.339 -0.389 1.090
## urg_medium 0.175 0.484 0.363 0.720 -0.819 1.170
## urg_high 0.077 0.397 0.194 0.848 -0.739 0.892
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Determine the potential moderating effect of urgency as continuous moderator. High urgency is chosen as inttercept category
urgency_cont_w8 <- rma.mv(effect_size, variance, mods = ~ urg, random = list(~1 | effectsize_id, ~ 1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes_filtered2)
summary(urgency_cont_w8, digits=3)##
## Multivariate Meta-Analysis Model (k = 29; method: REML)
##
## logLik Deviance AIC BIC AICc
## -29.872 59.743 67.743 72.776 69.648
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.464 0.681 29 no effectsize_id
## sigma^2.2 0.000 0.000 8 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 26) = 148.978, p-val < .001
##
## Test of Moderators (coefficients 2:3):
## F(df1 = 2, df2 = 26) = 0.068, p-val = 0.934
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.427 0.167 2.560 0.017 0.084 0.770 *
## urglow -0.077 0.397 -0.194 0.848 -0.892 0.739
## urgmedium 0.099 0.364 0.271 0.788 -0.649 0.846
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Determine the potential moderating effect of time budget to collision.
tbtc_w8 <- rma.mv(effect_size, variance, mods = ~ tbtc, random = list(~ 1 | effectsize_id, ~1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes_filtered2)
summary(tbtc_w8, digits=3)##
## Multivariate Meta-Analysis Model (k = 29; method: REML)
##
## logLik Deviance AIC BIC AICc
## -30.518 61.035 67.035 70.923 68.079
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.444 0.666 29 no effectsize_id
## sigma^2.2 0.000 0.000 8 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 27) = 150.829, p-val < .001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 27) = 0.028, p-val = 0.868
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.475 0.281 1.688 0.103 -0.102 1.052
## tbtc -0.005 0.028 -0.168 0.868 -0.061 0.052
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Further Moderator Analyses
Simulator fidelity
sim_fidelity_w8 <- rma.mv(effect_size, variance, mods = ~ sim, random = list(~1 | effectsize_id, ~ 1 | study_id), sigma2=c(NA,0), tdist=TRUE, data=effectsizes_filtered2)
summary(sim_fidelity_w8, digits=3)##
## Multivariate Meta-Analysis Model (k = 29; method: REML)
##
## logLik Deviance AIC BIC AICc
## -29.350 58.699 66.699 71.732 68.604
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.446 0.668 29 no effectsize_id
## sigma^2.2 0.000 0.000 8 yes study_id
##
## Test for Residual Heterogeneity:
## QE(df = 26) = 147.021, p-val < .001
##
## Test of Moderators (coefficients 2:3):
## F(df1 = 2, df2 = 26) = 0.567, p-val = 0.574
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.755 0.368 2.051 0.050 -0.002 1.512 .
## simlow fidelity -0.229 0.486 -0.472 0.641 -1.229 0.770
## simmedium fidelity -0.409 0.403 -1.016 0.319 -1.236 0.419
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#der_wo <- rma.mv(effect_size, standard_error, mods = ~ der, random = list(~1 | effectsize_id, ~ 1 | study_id), tdist=TRUE, data=effectsizes_filtered)
#summary(der_wo, digits=3)Interaction with other road users
iru_w8 <- rma.mv(effect_size, standard_error, mods = ~ iru, random = list(~1 | effectsize_id, ~ 1 | study_id), tdist=TRUE, data=effectsizes_filtered2)
summary(iru_w8, digits=3)##
## Multivariate Meta-Analysis Model (k = 29; method: REML)
##
## logLik Deviance AIC BIC AICc
## -29.828 59.657 67.657 72.840 69.475
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.195 0.442 29 no effectsize_id
## sigma^2.2 0.000 0.000 8 no study_id
##
## Test for Residual Heterogeneity:
## QE(df = 27) = 48.173, p-val = 0.007
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 27) = 0.193, p-val = 0.664
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.487 0.187 2.598 0.015 0.102 0.871 *
## iruyes -0.112 0.255 -0.439 0.664 -0.636 0.412
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Level of automated driving
lad_w8 <- rma.mv(effect_size, standard_error, mods = ~ lad, random = list(~1 | effectsize_id, ~ 1 | study_id), tdist=TRUE, data=effectsizes_filtered2)
summary(lad_w8, digits=3)##
## Multivariate Meta-Analysis Model (k = 29; method: REML)
##
## logLik Deviance AIC BIC AICc
## -29.831 59.662 67.662 72.845 69.480
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.196 0.443 29 no effectsize_id
## sigma^2.2 0.000 0.000 8 no study_id
##
## Test for Residual Heterogeneity:
## QE(df = 27) = 48.097, p-val = 0.007
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 27) = 0.263, p-val = 0.612
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.320 0.244 1.311 0.201 -0.181 0.820
## lad1 0.147 0.286 0.513 0.612 -0.440 0.734
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Publication Bias
Trim-and-fill method
##
## Estimated number of missing studies on the left side: 7 (SE = 3.5718)
##
## Random-Effects Model (k = 36; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -46.0221 92.0443 96.0443 99.1549 96.4193
##
## tau^2 (estimated amount of total heterogeneity): 0.6845 (SE = 0.1841)
## tau (square root of estimated tau^2 value): 0.8273
## I^2 (total heterogeneity / total variability): 90.23%
## H^2 (total variability / sampling variability): 10.23
##
## Test for Heterogeneity:
## Q(df = 35) = 259.2019, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.2021 0.1463 1.3809 0.1673 -0.0847 0.4889
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Egger Regression Symmetry Test
(eggers.test not working, bc ‘overall’ is no ‘meta’-type. Regtest is carried out with the two-level model)
##
## Regression Test for Funnel Plot Asymmetry
##
## model: mixed-effects meta-regression model
## predictor: standard error
##
## test for funnel plot asymmetry: z = 1.2719, p = 0.2034The p-value is > 0.05, publication bias can be rejected.
Another method for calculating Egger’s regression test - same result though.
egger <- rma.mv(effect_size, variance, mods = standard_error, random = list(~ 1 | study_id, ~ 1 | effectsize_id), tdist= TRUE, data=effectsizes_filtered2, method="REML")
summary(egger , digits=3)##
## Multivariate Meta-Analysis Model (k = 29; method: REML)
##
## logLik Deviance AIC BIC AICc
## -29.702 59.404 67.404 72.587 69.222
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.000 0.000 8 no study_id
## sigma^2.2 0.415 0.644 29 no effectsize_id
##
## Test for Residual Heterogeneity:
## QE(df = 27) = 145.893, p-val < .001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 27) = 1.618, p-val = 0.214
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt -0.358 0.635 -0.563 0.578 -1.660 0.945
## mods 2.839 2.232 1.272 0.214 -1.741 7.418
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Rank-Correlation Test
##
## Rank Correlation Test for Funnel Plot Asymmetry
##
## Kendall's tau = 0.1872, p = 0.1608